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The number-conserving quantum phase space description of the Bose-Hubbard model is discussed 
for the illustrative case of two and three modes, as well as the generalization of the two-mode case 
to an open quantum system. The phase-space description based on generalized SU(M) coherent 
states yields a Liouvillian flow in the macroscopic limit, which can be efficiently simulated using 
Monte Carlo methods even for large systems. We show that this description clearly goes beyond the 
common mean- field limit. In particular it resolves well-known problems where the common mean- 
field approach fails, like the description of dynamical instabilities and chaotic dynamics. Moreover, 
it provides a valuable tool for a semi-classical approximation of many interesting quantities, which 
depend on higher moments of the quantum state and are therefore not accessible within the common 
approach. As a prominent example, we analyse the depletion and heating of the condensate. A 
comparison to methods ignoring the fixed particle number shows that in this case artificial number 
fluctuations lead to ambiguities and large deviations even for quite simple examples. 

PACS numbers: 03.75.Lm, 03.65.-w 



I. INTRODUCTION 

The physics of ultracold atoms in optical lattices has 
made an enormous progress in the last decade, since it 
is an excellent model system for a variety of fields such 
as nonlinear dynamics or condensed matter physics. The 
dynamics of bosonic atoms can be described by the cel- 
ebrated Bose-Hubbard Hamiltonian, which is a paradig- 
matic model for the study of strongly correlated many- 
body quantum systems [l[ . Such systems are hard to deal 
with theoretically since the dimension of the respective 
Hilbert space increases exponentially both in the particle 
number and in the number of lattice sites. Therefore, 
approximations to the dynamics are of considerable in- 
terest. 

However, things can become surprisingly simple if the 
atoms undergo Bose-Einstein condensation. In many sit- 
uations, the mesoscopic dynamics of a Bose-Einstein con- 
densate (BEC) is extremely well described by the (dis- 
crete) Gross-Pitaevskii equation (GPE) for the macro- 
scopic wavefunction of the condensate (see, e.g., Q). 
This mean-field description is often referred to as 'classi- 
cal' since the macroscopic limit N — > oo of the dynamics 
in second quantization is formally equivalent to the limit 
h — > which yields classical mechanics as the limit of 
single particle quantum mechanics. In practice, this ap- 
proximation is frequently derived within a Bogoliubov 
approach, factorizing the expectation values of products 
of operators into the product of expectation values. As a 
consequence there is no obvious indicator to quantify the 
errors or to specify the scope of validity. Moreover, many 
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interesting observables which involve higher moments of 
quantum state are not accessible within this approxima- 
tion. A more detailed overview over established methods 
and a comparison to the approach presented here will be 
given in Sec. IVII1 

The phase space formulation of quantum mechanics is 
nearly as old as the theory itself and has a wide range of 
applications, especially in quantum optics. Beneath the 
illustrative insight into the dynamics provided by this 
formulation, many techniques and methods were devel- 
oped in this context. However, only a small fraction of 
the literature is dedicated to systems with intrinsic sym- 
metries like the Bose-Hubbard Hamiltonian. In this case 
the dynamical group is isomorphic to the special unitary 
group SU(M), with M being the number of lattice sites, 
reflecting the conservation of the total particle number. 
A general algorithm to construct phase space distribu- 
tion functions for systems whose dynamical group has 
the structure of an arbitrary Lie group such as SU(M) 
has been developed only nine years ago 0|. It is based 
on the concept of generalized coherent states introduced 
by Gilmore [H and Perelomov Starting from these 
states we have surveyed the mathematical foundations 
of the number conserving phase space description and 
derived the exact evolution equations for the Husimi Q- 
function and the Glauber-Sudarshan P-function of the 
M-site Bose-Hubbard model in a preceding paper Q. 
One important consequence of the use of SU (M) coher- 
ent states is the different topology of the phase space, 
which is now isomorphic to the 2M — 2 dimensional Bloch 
sphere and therefore compact. Moreover, the phase space 
description based on SU(M) coherent states allows for 
a deeper analysis of the many-particle-mean-field corre- 
spondence, since these states are of high physical signif- 
icance. As they are equivalent to the fully condensed 
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product states, they provide an excellent tool to describe 
and study derivations from the macroscopic state which 
also determine the region validity of the mean-field ap- 
proximation. 

In this paper we will illustrate the methods originally 
developed for the M-site system Q for the instructive 
case of small Bose-Hubbard systems. Therefore we con- 
sider a two-mode Bose Hubbard model 

H = - /S.(t) (^a\a 2 + a\a^ + e 2 a\a 2 + t\a\ai 

+ j(a\ 2 al + a 2 2 at), (1) 

and its generalization which takes into account elastic 
collisions with background gas leading to phase noise de- 
scribed by a master equation. Such a system can be real- 
ized experimentally by confining a BEC in a double-well 
trap @, H, H, E3|. Moreover, we also consider the case of 
an explicitly time-dependent tunneling-element A, which 
leads to mixed regular-chaotic mean- field dynamics. This 
can be easily implemented in an experimental setup via 
a variation of the intensity of the laser forming the opti- 
cal potential. However, this model can also describe the 
fundamental phenomena of two weakly coupled BECs in 
more general setups and Landau- Zener transitions 
between different Bloch bands [H, EH . Early theoreti- 
cal studies of the system dynamics without an external 
influence and a time-dependant driving were reported in 
[OIIIEI, EH- Furthermore, we use the methods presented 
here to analyse the interesting and yet not completely un- 
derstood case M = 3, i.e. the three-mode Bose-Hubbard 
system: 

H = — A12 (a\a 2 + a^aij - A23 (a5,a 3 + a\a 2 ^j 

3 3 

i=i 3=1 

The realization of such a linear triple- well trap for a BEC 
is not as straightforward as in the case of a double-well 
trap, but nevertheless possible in optical setups [lTj or on 
an atom-chip. Besides it is a prominent model system be- 
cause the mean-field dynamics is classically chaotic (see, 

e.g. miiaiEiij). 

In particular, this paper is organised as follows: In the 
sections [TTI and Hill we briefly recall the main results from 
the number conserving phase space approach and pro- 
vide some illustrative examples. Then follows an analysis 
of the exact dynamics of the quasi probability distribu- 
tions in Sec. HVl Here, we are especially interested in 
the macroscopic limit, since the evolution equations con- 
verge to a classical Liouville equation for large particle 
numbers. This result enables us to go beyond the usual 
mean-field limit which considers only point-distributions 
and to enhance the region of applicability of the approx- 
imation to arbitrary initial states. To emphasize this, we 
study in Sec.|V]the behaviour around dynamical instabil- 
ities where the common mean-field approach is known to 



fail badly. In the case of an isolated instability this effect 
is known as the breakdown of mean- field. Here, we show 
that the Liouville dynamics not only resolves this break- 
down due to an isolated instability, but also provides an 
ideal tool to describe systems where the mean-field dy- 
namics faces large chaotic regions, as for example in a 
driven two mode system. 

Having demonstrated the considerable advantages of 
the Liouville description, we are now able to address novel 
questions, like the dynamics of the three-mode system 
in Sec. I VII While being conceptually quite simple and 
much easier to implement than other methods (cf. e.g. 
[H, HH) the Liouville description provides a valuable tool 
to estimate the depletion of the condensate mode and to 
infer the nature of the many-body quantum state. The 
comparison to exact numerical results show clear agree- 
ment with the results of the Liouville equation while the 
common mean-field approach fails. 

Since the methods developed in Q are not restricted 
to the Bose-Hubbard model, but can be applied to ev- 
ery problem within the same symmetry group, there is 
a large variety of possible applications. As one exam- 
ple, we discuss the heating of a two-mode BEC result- 
ing from elastic scattering with the background gas in 
Sec. lVI CI While our methods predict the proper results, 
namely a quasi-classical set of Langevin equations in the 
mean-field approximation, they allow to go even further 
and predict many interesting effects such as the decay 
of the coherence factor and the increase of the variances, 
which are not accessible within the usual mean-field limit. 
This again underlines the benefit from the Liouville de- 
scription and the distinction to common mean-field ap- 
proaches. 

In Sec. lVIIl we then give a short overview over different 
approaches to derive and extend a mean-field description 
for the systems studied here. Especially, we compare our 
approach based on the dynamical symmetry and there- 
fore explicitly including the conservation of the parti- 
cle number to common phase space methods (cf. e.g. 
[HIHH). One important result is that calculations based 
on a U(l) symmetry breaking description using common 
Glauber coherent states show deviations from the exact 
dynamics which are significantly larger and appear on a 
shorter time scale compared to the number conserving 
approach. 



II. ALGEBRAIC STRUCTURE AND 
GENERALIZED COHERENT STATES 

The number conserving phase space description of the 
Bose-Hubbard dynamics has been introduced in Q , start- 
ing from the generalized coherent states of Gilmore [4] 
and Perelomov [f| . Here we will only recall the main re- 
sults for the special case of two and three modes for the 
sake of completeness. 

Generalized coherent states for arbitrary dynamical 
Lie groups are defined by the action of a translation 
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operator onto a reference state. For instance, the cel- 
ebrated Glauber coherent states for the Heisenberg-Weyl 
group H4 are defined by the relation \a) — D(a)\Q) = 
£ aa 1 -a °|o) ; where the reference state is just the vacuum 
state. The parameter space of the translation operators 
D(a) and thus the space of the coherent states is isomor- 
phic to the classical phase space C ~ K x 1. 



A. The two-mode case 

The situation is different for a BEC in a double well 
trap since the dynamical group is now SU (2) spanned by 
the angular momentum operators 



J X 


1 

- [a[a 2 + 




Jy ~ 






Jz = 


- (a\a 2 - 


a\ai 



(3) 



The Hamiltonian {!]) then can be rewritten as 

H = -2AJ X + 2eJ z + UJ 2 Z (4) 



with e = € 2 — (\ up to a constant term [lj, [15|, l26l l55l|. 
The Casimir operator J 2 = N/2 (N/2 + 1) reveals the 
relation between the algebraic structure and the total 
particle number N — a\ai + a\a,2- 

Generalization of the concept of the translation opera- 
tor to the su(2) algebra yields an rotation operator, such 
that the SU(2) or so called Bloch coherent states are 
defined as follows: 

\6,<t>) = R{6,<t>)\N,Q) (5) 

_ g—i9[Jx sin cp—Jy cos cp) I jy 
l 

= E ( N Y co S (i) ni sm(i) n2 e- i ^1'\m,n 2 ). 

ni +n 2 =N ^ n2 ' / 

The group of rotations R(9, <fi) and therefore the param- 
eter space of coherent states can be described by two 
angles 0<6<tt,0<(P<2tt and is thus isomorphic to 
a Bloch sphere S 2 . In physical terms, this can be under- 
stood as follows: The z-component of the Bloch vector 
describes the population imbalance of the two modes, 
while the polar angle represents the relative phase. Since 
this variable is cyclic and not defined if all particles are 
in a single well (i.e. at the poles), the topology is clearly 
that of a sphere. 

Note that the SU(M) coherent states are equivalent to 
the product states, i.e. they represent a pure condensate 
[USB For two modes, M — 2, this relation simply reads 
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cos(f)al +sin(f)e-^4) |0,0). (6) 



JV 



is equal to the relative phase between the modes and 
p = cos 2 (9/2) denotes the population of the second mode. 



B. The three- mode case 

For the three-mode system the dynamical group is iso- 
morphic to SU(3) which is spanned by the eight genera- 
tors 

X\ = a\a,i — a\a,2, 



1 ' a\di + a\a,2 — 2030,3 



X 2 = 



Yt, = / ! tt\.a . — a' rj a k 

Zk = a' k aj + aja fe , (7) 

with k — 1,2,3 and j = k + 1 mod 3, which form a 
generalized angular momentum algebra (cf. e.g. [(j [28l|). 

A convenient choice of basis is given by the eigenstates 
of the Casimiroperator 4N(N/3 + l) with N — Yli=i 4^* 
and the two elements of the Cartan subalgebra Aj and 
X 2 with 

hk\n 1 ,n 2 ,n 3 = N - n x - n 2 ) = n k \ni,n 2 ,n s ). (8) 

In this basis the reference state is chosen to be the highest 
weight state |iV, 0,0) and the generalized coherent states 
I ft) are obtained by the natural extension of (0 to higher 
dimensional rotations. Note, that SU(3) is a two instead 
of one parameter Lie group, thus the parameter space 
of generalized the coherent states is now isomorphic to 
SU(3)/U(2) and therefore four-dimensional. In order to 
visualize the phase space dynamics we must hence use 
projections and Poincare sections, which are not neces- 
sarily be unique. For further mathematical details, see 
e.g. [Ill for SU(3), as well as |g,|25| for the general case 
SU(M). 

However, for any actual calculation or numerical sim- 
ulation the equivalent parametrization of the generalized 
coherent states |ft) via the complex mode amplitudes ipi 



|ft) = 




0,0,0) 



(9) 



This representation of SU(2) coherent states reveals a 
parametrization in more physical terms p, q, where q = 4> 



proves to be more practical. Using the conservation of the 
total particle number which corresponds to the normal- 
ization 2 = |V'i| 2 + |V'2| 2 + |V'3| 2 = 1 , we can fully char- 
acterize the state vector by the occupations P2 = \ip 2 \ 2 , 
Ps = IV'sl 2 an d the relative phases q x = arg(?/>2) — arg(V'i) 
and qs — arg(^2) — argf^), as discussed in the case of the 
two mode system. Therefore we will use this parametriza- 
tion in the following. 

Note, that such a description in terms of SU(3) co- 
herent states has already proven useful in the case of 
a circular arrangement with uniform tunneling elements 
Aij = A = const. [3(| and in the study of quantum 
discrete vortices [3l[. A comparison of these results to 
calculations based on common Glauber coherent states 
can be found in [32j |. 




FIG. 1: Husimi functions in a Mercator projection of (a) the 
coherent state \p = 1/2, q = 0), (b) the Fock state \m — 
N/2, ri2 = N/2) and (c) the ground state of the Hamiltonian 
|T]) for A = 1, U — 10 and e = 0. The number of particles is 
N = 40 in all plots. 




FIG. 2: (Color online) Husimi density of the cat states 
\ip±) oc |0.2, 0) ± 0.8, 0) (a),(b) and the incoherent sum 
p+ = ([0.2,0) (0.2, 0| + |0.8,0) (0.8, 0|)/2 (c) for N = 20 par- 
ticles. The zeros of the Husimi function are marked by red 
crosses. 



III. PHASE SPACE DISTRIBUTIONS 

With the help of the generalized coherent states, which 
will be denoted |0) in the following (independently of 
the corresponding dynamical group) , one can readily in- 
troduce quasi phase space distributions for a BEC in a 
double or triple well trap. The Glauber-Sudarshan P- 
distribution is defined as the diagonal representation of 
the density operator p in generalized coherent states 

p = I p(Q)\n){n\dfj,(n), (io) 

where 51 represents the respective parametrization and 
dp(Q) stands for the invariant measure on the respective 
phase space, denoted by X. Due to the overcompleteness 
of the coherent states, the P-function does always exist 
but is usually not unique. Furthermore it is not positive 
definite and often highly singular. On the other hand, 
the Husimi Q-function defined as the expectation value 
of the density operator in generalized coherent states, 

Q(Q) = (Sl\p\Q), (11) 

is unique, regular and positive definite. Thus the Q- 
function is especially suited for illustrations, while both 
quasi distribution functions will be used for actual cal- 
culations. However, the Q-function is also not a prob- 
ability distribution function in the strict sense since it 
does not give the correct marginal distributions. Note 
that it is also possible to define the Wigner function on a 
spherical phase space but this leads to much more com- 
plicated expressions than the P- and Q-function, since 
its construction uses harmonic functions on the respec- 
tive phase space. This makes actual calculations a hard 



task even for two lattice sites corresponding to SU(2) 
(see, e.g., the contradictory results in [33[ and (34|) and 
almost impossible for larger systems. Furthermore, a pos- 
itive P-representation respecting SU(M) symmetry has 
been introduced and analyzed in [35| . 

To get a first impression of the phase space distribu- 
tions considered here we have plotted the Husimi func- 
tion of the two mode system for (a) the coherent state 
\p = 1/2, q = 0), (b) the Fock state |ni = N/2, n 2 = N/2) 
and (c) the ground state of the Hamiltonian (JXJ) for A = 1 
and U — 10 in a Mercator projection of the sphere in 
Fig. [TJ The coherent state is maximally localized at the 
position [p, q) and thus closest to a point in classical 
phase space. On the contrary a Fock state is localized 
around p — rii/ 'N . The phase q is completely delocal- 
ized in the sense that the Husimi function is uniform in 
q. Finally one can clearly visualize the number squeez- 
ing of the eigenstate (c) in comparison with the coherent 
state (a) due to the interaction term in the Hamiltonian 
- a fact which is desirable for matter wave interferometer 
experiments [36j]. 

The overcompleteness of the basis of (generalized) co- 
herent states guarantees that one can uniquely recon- 
struct the quantum state out of the quasi probability 
distribution. It is even possible to determine a pure 
state solely from the N zeros of the Husimi function 
or the Bargmann function up to a global phase factor 
[37l ]. Thus these zeros carry the essential information 
about the quantum state and its coherence properties. 
As an example, Fig. [2] compares the Husimi distribution 
of a cat state in the two mode case, i.e. a coherent su- 
perposition of two SU{2) coherent states \ip±) oc \p — 
0.2, q = 0) ± \p = 0.8, 9 = 0), with the incoherent sum of 
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these two states p+ = (|0.2,0) (0.2, 0| + 10.8,0) (0.8,0|)/2. 
The global shape of the Husimi function of the three 
states appears quite similar with significant differences 
only around the point (p, q) — (0.5,0). The Husimi den- 
sity is increased for the coherent superpositions |V-'+) m 
(a) and decreased for \ip-} in (b) in comparison with the 
incoherent sum, reflecting constructive or destructive in- 
terference, respectively. The complete information about 
the pure states \ip±) is coded in the distribution of the 
N = 20 zeros of the Husimi function, which are plotted 
as red crosses in the figure. The Husimi function of the 
mixed state /3+ has no zero at all, whereas the Husimi 
function of a coherent state (compare also Fig. [1]) has a 
single iV-fold degenerate zero opposite to the maximum 
of the quasi probability distribution. 




1 2 1 2 



q/rc q/n 

FIG. 3: (Color online) Classical phase space structure: Lines 
of constant energy Ti(p,q) for A = 1, e = 0, and g — (left) 
resp. g — 10 (right). 



IV. DYNAMICS AND THE MACROSCOPIC 
LIMIT 



For the P-function we need an associated 2? € -algebra de- 
fined by an integration by parts: 



Evolution equations for the classical phase space dis- 
tributions can be derived using the mapping of operators 
in Hilbert space to differential operators on the classical 
phase space. The 2? £ -algebra representation of an arbi- 
trary hermitian operator A is defined such that the result- 
ing differential operator acting on the parametrization of 
the generalized coherent state projector has exactly the 
same effect as the original operator: 



A\Q) (n\ = v l (A)\n) , 
\n) (n\A = v e (Ay\n) (o| . 



(12) 
(13) 



Here, the superscript £ refers to side (left) from which 
the operator acts. The P^-algebra representation is well 
known for the Heisenberg-Weyl algebra. It has been in- 
troduced for the su(M)-algebra by Gilmore and cowork- 
ers (see [4] and references therein). Here, we just state 
the expression for the generators of the su(2)-algebra in 
the angular parametrization ([5]): 



{P(n)V e (H)}\n) (0|d/x(0) 



A' 



{f> e (H)p(n)}\n) (fi| <f/x(fi). (16) 



A 



Using this algebra, the evolution equation for the P- 
function can be derived via 



,, = i ^P(n)|n)(n|d/i(n) 



(17) 



v 1 {h) - v 1 {h)* p(n)\n) (n\ dn(n). 



A 



The T> -algebra representation of the relevant oper- 
ators and the evolution equation for the M-site Bosc- 
Hubbard-Hamiltonian have been calculated in the pre- 
ceding paper [f|. They are generally applicable to every 
problem within this symmetry class. 



J<l> 



N 



■ sm f 



N 



i d 
2 COt 2^ H 

1 9 d 
tan — — — 

2 2 86 



cos 



9 d 
289 



2 89 



N 



1 • n d 
2 Sm9 d9 



i 8 
280 



(14) 



Independently of the number of modes M, the evo- 
lution equation for the Q-function can be expressed in 
terms of the differential operators T> e : 



tr(/5|Q) (Q|) 

-i tr(pH\n) (Q| - \ (Q\ Hp) (15) 
-i(V e (H)-V e (Hy)Q(Q). 



A. Mean-field dynamics of the two-mode system 

Since the mapping to the differential operators is lin- 
ear, the explicit evolution equations for the Q- and P- 
function can be directly calculated from the imaginary 
part (cf. Eqn. (fT5"|) and (jTTJ)) of the T> e operators (see 
Eqn. p^)l ). This yields the following results: 



8_ 

di 



Q(p,q) 



2e- 



8 



(18) 



+A 2y / p-p 2 sin 



8 1 - 2p 8 \ 

\ d P VP~ P dq J 

-u(n(1-2 P )-2 P (1-p)^) 



FIG. 4: Husimi density of the eigenstates \E n ) with n = 
1,15,23,34 (a-d) of the Hamiltonian © for A = 1, e = 0, 
UN — and N = 40 particles. Dark colors encode high 
values of the Husimi density. 



for the Husimi Q-function and 



2 — 

dq 



(19) 



+A 



p 2 sin q 



0_ 

'dp 

+U ({N + 2)(l-2p) + 



l-2p 



■. cos q- 



Vp-p 2 



P(p,q) 



for the Glauber-Sudarshan P-function. 

One can directly show that both evolution equations 
conserve the normalization. Furthermore it can be shown 
that the interaction terms ~ U cancel exactly for N = 1, 
since a single Boson obviously does not interact. The 
expectation values of the angular momentum operators 
([3]) are exactly given by the phase space averages [6| 



(J fc ) = {N + 2) / s k (p, 1)Q(P, Q)MP> l) and 



(J fe ) = N s k (p,q)P{p,q)dp{p,q), 
Js 2 

where k = x,y,z and 

/ y/p(l -p) cos(g) 
S= y/p(l-p) sm(q) 



(20) 



(21) 



is the classical Bloch vector. 

These exact equations provide a valuable starting point 
to study the limit of large particle numbers, since it 
emerges quite naturally in the phase space picture. In 
the macroscopic limit N — > oo with g = UN fixed, the 
second order derivative terms in the evolution equations 
become negligible, as they are vanishing as 0(1/N). The 
appearance of a non-positive diffusion term as the differ- 
ence between the exact quantum dynamics and the classi- 
cal dynamics has already been discussed in (38j . Thus, in 



FIG. 5: As Fig. H however for UN = 10. 



the macroscopic case one obtains a Liouville equation for 
a quasi-classical phase space distribution function p(p, q): 



dH d _ dH d 
dp dq dq dp 

Wp, q),p(p, q)} 



pip, q) 



(22) 



with the Poisson bracket {•,•} and the Gross-Pitaevskii 
Hamiltonian function 



H 



-2ep - 2Ay/p{l-p) cos(q) + f (1 - 2pf . (23) 



The macroscopic interaction strength depends upon the 
operator ordering and is thus given by g = UN if we 
start from the Q-function and by g = U (N + 2) for the 
P- function. This difference also vanishes in the macro- 
scopic limit N — > oo with g = UN fixed. Therefore the 
phase space picture gives rise to a Liouvillian flow for the 
time-evolution of the quasi probability distribution in the 
macroscopic limit. 

The common mean- field limit is subject to even more 
severe restrictions. To circumstantiate this point, we will 
have a closer look at the relation between the Liouville 
dynamics and the GPE. To obtain the latter, one must 
not only take the macroscopic limit of the evolution equa- 
tions (fT5)) and ([TU|) . but also assume a pure condensate, 
respectively a SU(M) coherent state which is maximally 
localized during the whole time evolution. In this case 
the dynamics of a phase the space distribution of the 
initial state is approximated by a point distribution, re- 
spective a single trajectory instead of an extended dis- 
tribution. Therefore, any information about quantities 
involving higher moments of the state is lost. 

So, if we assume a delta-distribution p(p' , q') = 5{p' — 
p)5(q' — q), we can derive the GPE starting from the Li- 
ouville equation (|22p . This yields the canonical equations 
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of motion 



P 



an 



-2Ay/p(l-p) sm(q) 



(24) 



q = ^ = -2e - A LJ^_ cos(g) - 3 (1 - 2p). 
dp Vp(1-p) 

These equations are equivalent to the celebrated discrete 
Gross-Pitaevskii equation 



d f tp 



'dt I tp: 



-A 



-A -< i ^' (25) 



if one identifies ipi = — p and ^2 = y^pexp (— iq) and 
neglects the global phase. In the non-interacting case 
[7 = 0, these mean-field equations of motion are exact 
and an initially coherent state remains coherent in time. 
Hence, no information is lost, since the higher moments 
can be easily reconstructed. However, in the interact- 
ing case an initially pure BEC typically deviates quite 
rapidly from the set of maximally localized states. More- 
over, by approximating the state by a product state one 
looses the entire information about the higher moments 
of the state. Therefore many quantities of interest as for 
example variances etc. are not accessible, whereas they 
can be easily estimated using the Liouville dynamics and 
phase space averages similar to Eqn. (|2"0|) . 

The classical phase space structure induced by the 
Hamiltonian function ([23]) is shown in Fig. [3J for A = 1 , 
e = and two different value of the interaction strength, 
<7 = and g = 10. In the linear case one recovers the 
simple Rabi or Josephson oscillations. One of the elliptic 
fixed points bifurcates if the interaction strength exceeds 
the critical value g = 2A, leading to the self-trapping 
effect 0EEE3- 

Furthermore, we have plotted the Husimi function for 
some selected eigenstates of the Bose-Hubbard Hamilto- 
nian fTj for the linear case UN = in Fig. 2] and for 
UN = 10 in Fig. \5\ The Husimi functions of the eigen- 
states localize on the classical phase space trajectories of 
the respective energy as shown in Fig. [3J where their sum 
covers the classical phase space uniformly, 



^2Q\ n )(p,q) = i- 



(26) 



In the linear case UN = 0, the classical trajectories cor- 
respond to the Josephson oscillations between the two 
wells as shown on the left-hand side of Fig. [3] This be- 
havior is directly mirrored by the eigenstates as shown 
in Fig. 31 The lowest eigenstate localizes on the classical 
fixed point (p, q) = (0, 0) while the others correspond to 
Josephson oscillations. 

The situation is more interesting, when the interaction 
strength exceeds the critical value for the self-trapping 
bifurcation, as shown exemplarily for UN = 10 in Fig. [5] 
For low energies one still finds Rabi modes, while the 
eigenstates of higher energies correspond to self-trapping 
trajectories with a persistent particle number difference 



and a running phase (cf. Fig. [3J right). However, the 
eigenstates always localize equally on the trajectories 
with p > 1/2 and p < 1/2 because of the symmetry of 
the Hamiltonian Most interestingly, quantum states 
can also localize at the unstable hyperbolic fixed point as 
shown in Fig. [5] (c). 



B. Mean-field dynamics of the three-mode system 

As in the case of the two-mode system, the exact evolu- 
tion equations of the phase space distribution functions 
can be calculated from the V 1 operators. The results 
given in [fj are considerably more complicated as in the 
case of only two modes. However, in the macroscopic 
limit one again recovers a Liouville equation 



d_ 

at 



-{H,p(pi,p 3 , <7i,93)} (27) 



with the macroscopic Hamiltonian function 



H 



eiPi + e 3 P3 + ^(Pi+pI + ( 1 -Pi- P3) 2 ) 
-2Aia\/i - pi - p 3y /p~[ cos(qi) 
-2A 23 ^1 -p! - P3\/p3COs(q 3 ). 



(28) 



The main new feature of the dynamics in three-levels 
systems is that the mean- field dynamics can become clas- 
sically chaotic (cf. e.g. [H, [13, HU and references 
therein), where the transition goes from regular dynam- 
ics for g = 0, to mixed chaotic dynamics for stronger 
nonlinearities \g\, whereas the mean- field dynamics of the 
two mode system with constant parameters is always reg- 
ular. Figure [HI (a) shows a Poincare section of the cor- 
responding phase space for g = — 5 and 7i = 1, where 
the dynamics is still mostly regular. However, a classi- 
cally chaotic strip is found between the elliptic fixed point 
a t (^3,93) = (0.66,0) and the running phase modes be- 
low. Chains of large regular islands are found within 
this chaotic sea (colored in green in the figure). An- 
other example is shown in Fig. [6] (b) for g = —10 and 
H = — 1. Again we observe elliptic fixed points with large 
values of p 3 separated from the running phase modes by 
a chaotic region in phase space. However, the regular 
running phase modes now occupy most of the energet- 
ically allowed region of phase space. Furthermore one 
observes a novel fixed point embedded in this region of 
phase space. 

This phase space structure is typical for the route to 
chaos in Hamiltonian systems described by the KAM and 
the Poincare-Birkhoff theorem. The dynamics becomes 
chaotic in the phase space region close to the separatrix. 
Resonant tori are destroyed and decompose into a series 
of elliptic and hyperbolic fixed points, where the elliptic 
ones form the island chains embedded in the chaotic re- 
gion. The remaining invariant tori confine this chaotic 
sea. 
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FIG. 6: (Color online) Poincare section of the classical phase 
space (gi = 0, see text for details) for ej = 2, 63 = 4, A12 = 
A23 = 1, g = — 5 and a fixed total energy 7i = +1 (a) and 
—g = 10 and Ti = — 1 (b). Trajectories around the elliptic 
fixed points are coloured in red to guide the eye. 



Note that the use of a Poincare section is not straight- 
forward for the dynamical system considered here. As 
usual, the variables (^3,93) are plotted when the trajec- 
tories intersect the q\ = 0-plane with qi > 0. The re- 
maining dynamical variable p± follows from energy con- 
servation. However the uniqueness has to be checked ex- 
plicitly. For example, the assignment of p\ is not unique 
in Fig.[6](b) in a small area of phase space which is shaded 
in grey in the figure. Solving equation [)28p for p\ gives 
two possible solutions. In order to make the calculation 
of pi unique we define that only the solution with smaller 
value of pi is plotted in the Poincare section. Further- 
more, large regions of the (j>3, 53) -plane are energetically 
forbidden for qi = (cf. Fig. [5]). 

Let us have a closer look at the mixed classically- 
chaotic phase space shown in Fig. [5] (b). The dynam- 
ics is also shown in a three-dimensional phase space plot 
in Fig. [7J Only few trajectories are plotted exemplar- 
ily. The periodic orbits (coloured in red) correspond to 
the elliptic fixed points in the Poincare section (Fig. [5]) . 
Two periodic orbits are found on the left-hand side of 
the figure, oscillating with a small amplitude both in the 
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FIG. 7: (Color online) Upper panel: mean-field phase space 
trajectories for ei = 2, £3 = 4, A12 = A23 = 1, g = —10 
and a fixed total energy H = — 1. Lower panel: Time series 
(P3{t) = \ip3 (t)\ 2 ) and power spectrum of a trajectory in the 
chaotic sea. 



populations and in the relative phases. Another periodic 
trajectory shows a larger amplitude of the oscillations of 
the populations pi and ps and a running phase 53 . This 
periodic trajectory is embedded in a set of similar quasi- 
periodic ones, of whom one is plotted as a blue line in 
the figure. Furthermore, a trajectory in the chaotic sea 
is plotted in grey. It is observed that the chaotic sea splits 
into two distinct regions with narrow links: One around 
the periodic orbits on the left and one in the vicinity of 
the running-phase trajectory on the right-hand side. 

This 'bottleneck' effect leads to an intermittent dy- 
namics: A trajectory can be confined to certain regions 
in the chaotic sea for a long time, where its time evolution 
looks quasi regular. Chaotic bursts to the remaining part 
of the chaotic sea are found between the quasi-regular 
oscillations. This is shown in Fig. [7J where the time se- 
ries P3(t) for one particular trajectory in the chaotic sea 
is plotted as well as the corresponding power spectrum. 
One observes long quasi-regular oscillations with a small 
amplitude intermittent by large amplitude bursts. Phys- 
ically this means that the BEC mostly stays in the right 
well for a long time with irregular bursts to the middle 
well. 

In conclusion, we stress that in any case, two steps 
of approximation were necessary to derive the Gross- 
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FIG. 8: Breakdown of the mean-field approximation around a hyperbolic fixed point in phase space. Shown is the dynamics 
of an initially coherent state located at po = 0.9045 and go = at times t = 0,0.12,0.24,0.36,0.48 (from left to right) for 
A = 1, UN — 10 and N — 20 particles. The exact quantum evolution of the Husimi function (top) is compared to the classical 
Liouvillian dynamics (middle) and to the dynamics of a classical phase space ensemble (bottom). 



Pitaevskii equation: In the first step we have neglected 
the quantum noise term in the evolution equations (|18[) 
and (|19p which is always possible in the macroscopic limit 
N — > oo with g = UN fixed. Only in the second step to 
derive the GPE one assumes a nearly pure condensate, 
so that its phase space representation is strongly local- 
ized and can be approximated by a point distribution or 
a single trajectory in phase space. 

Thus it is not only possible to simulate the dynamics of 
a strongly localized product state, but also of every other 
possible initial state using the Liouville equation while 
the error vanishes as \/N . Whereas the single-trajectory 
mean- field approximation will break down if the quantum 
state differs significantly from a product state, the phase 
space description still gives excellent results. This issue 
and its implications will be further investigated in the 
following. 



V. ANALYSIS OF DYNAMICAL 
INSTABILITIES AND CHAOS 

The phase-space description of the Bose-Hubbard 
model is especially suited to explore the correspon- 
dence of the quasi-classical mean-field approximation and 
the many-particle quantum dynamics. The established 
derivations of the Gross-Pitaevskii equation assume a 
fully coherent quantum state, so that the mean-field ap- 
proximation is restricted to this class of quantum states. 
In fact it has been shown that the approximation is no 
longer valid if the Gross-Pitaevskii dynamics becomes 
classically unstable [12, [Hj]. In this section, we inves- 



tigate the behaviour of the mean-field approximation 
around localized, as well as extended regions of dynam- 
ical instabilities. It is shown, that the Liouville dynam- 
ics not only resolves well-known problems of the com- 
mon mean-field description, but enables us to even study 
highly chaotic systems. 



A. Resolution of the break-down of mean-field 

A particular illustrative example of the failure of the 
common mean-field equations was introduced by Anglin 
and Vardi [26|, [39j], who demonstrated the breakdown 
of the mean-field approximation for a two-mode BEC 
around the hyperbolic fixed point shown in Fig. [3] on the 
right-hand side. Therefore, this is an excellent example 
to demonstrate that the approximation by the Liouville 
dynamics clearly goes beyond the usual mean-field de- 
scription and to enlight the reason of the break-down. 
The top row of Fig. [5] shows the resulting evolution of 
the quantum Husimi Q-function. Initially the system is 
in a fully condensed state, i.e. a SU(2) coherent state 
at po — 0.9045 and go = 0, and thus the Husimi func- 
tion is maximally localized. In the course of time the 
condensate approaches the hyperbolic fixed point, where 
the Husimi function rapidly diffuses along the unstable 
classical manifold. The coherence is lost and the Gross- 
Pitaevskii equation is no longer valid. This is further 
illustrated in Fig. [9l where the evolution of the quantum 
Bloch vector (J)/N (solid blue line) is compared to its 
classical counterpart s (solid red line). Both agree well in 
the beginning, when the system is fully condensed. How- 
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FIG. 9: (Color online) Time evolution of the Bloch vec- 
tor (J) /N approaching the hyperbolic fixed point (solid blue 
line) for the same parameters as in Fig. [8] The ensemble 
average over 500 trajectories (dashed green line) closely fol- 
lows the quantum result, while a single trajectory (solid red 
line) breaks away in the vicinity of the hyperbolic fixed point 
(marked by a cross). 



ever, the mean-field approximation breaks down as they 
approach the hyperbolic fixed point (marked by a cross 
in the figure). The quantum Bloch vector penetrates into 
the Bloch sphere, while the classical vector s is restricted 
to the surface. 

However, this breakdown of the mean-field approxi- 
mation is only due to the neglect of higher moments of 
the quantum state and not a consequence of a failure 
of the quasi-classical Gross-Pitaevskii dynamics. Thus 
it is easily resolved in quantum phase space. The mid- 
dle row of Fig. [5J shows the evolution of a classical phase 
space distribution propagated by the Liouville equation 
(|2"2")) which coincides with the Husimi function at t = 0. 
With increasing particle number N the width of the ini- 
tial distribution decreases such that one recovers a sin- 
gle phase space point in the macroscopic limit N — > oo. 
One observes that the classical distribution captures the 
essential features of the quantum evolution - the spread- 
ing along the unstable manifold of the hyperbolic fixed 
point. Due to the neglect of the second order differen- 
tial term in Eqn. (fT5|) this spreading is a little overes- 
timated while the quantum spreading in the orthogonal 
direction is absent. Thus, this can be construed as a 
disregard of quantum noise which guarantees the uncer- 
tainty relation. Furthermore we can interpret the Liou- 
villian flow in terms of a classical phase space ensemble 
as shown in the bottom row of Fig. [51 At t = this en- 
semble is generated by 500 phase space points to mimic 
the quantum Husimi distribution. Afterwards all trajec- 
tories evolve according to the Gross-Pitaevskii equation 
(I25p resp. (124|) . This approximation is comparable to nu- 
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FIG. 10: Stroboscopic Poincare section of the classical dy- 
namics of a driven two-mode BEC for Ao = 1, g = 5, u = 2tt 
and Ai = 0,0.2, 0.5 (from top left to bottom). 

merical Monte Carlo methods and especially suited for 
the generalization to larger systems, since it is easily im- 
plcmentablc. Fig.[9]shows the quasi-classical expectation 
value of the Bloch vector (dashed green line) calculated 
by the simple phase space average (|20[) in comparison to 
the quantum Bloch vector (3)/N. One observes an ex- 
cellent agreement. Note that the expectation values cal- 
culated from the classical distribution function and the 
phase space ensemble are indistinguishable on this scale 
and thus only one curve is shown in the figure. 

B. Description of chaos in a driven two- mode 
system 

Periodically driven Bose-Hubbard systems are of cur- 
rent interest (see, e.g., [4(| and references therein). In 
this paragraph, we consider the dynamics in a double- 
well trap for a time-dependent tunneling matrix element 

A(t) = A + Ai cos(wi). (29) 

Such a driving can be implemented easily in an opti- 
cal setup by varying the intensity of the laser beams 
forming the optical lattice. Figure [TU] shows a strobo- 
scopic Poincare section of the mean-field dynamics for 
the parameters Aq = 1, g = 5, Ai = 0,0.2,0.5 and 
oj = 2ir so that the period is T = 1. The fixed points 
of the stroboscopic mapping correspond to the nonlin- 
ear Floquet solutions of the driven GPE (cf. also [40j| ) . 
For Ai = one obviously recovers the undriven system, 
which is completely regular. A mixed regular-chaotic dy- 
namics is found when Ai is increased, the region around 
the hyperbolic fixed point (p,q) = (1/2, n) being the 
first to become chaotic. Fig. [10] (b) shows the phase 
space for Ai = 0.5, where only tiny regular islands 
have survived, among them regions around the symmet- 
ric fixed point (p, q) = (1/2, 0) and the self-trapping fixed 
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FIG. 11: Comparison of quantum and classical dynamics in phase space in the chaotic sea for UN = 5, Ao = 1, u) = 2tt and 
Ai = 0.5 and N = 20 particles. Upper panels: Exact Evolution of the Husimi density at times t = 0,1,2,3,4 (from left to 
right). Lower panels: Dynamics of an ensemble of 150 classical trajectories. 



points. The regular islands around (p, q) = (l/2,7r/2) 
and (p, q) — (l/2,37r/2) correspond to a period doubled 
fixed point. They are of particular interest since they 
support regular oscillations of the population with an ex- 
tremely large amplitude as illustrated in Fig. [TJJ 

If the classical dynamics is unstable, an initially co- 
herent quantum state will be strongly distorted. For an 
isolated unstable fixed point, the state spreads along the 
unstable manifold while it is squeezed in the orthogonal 
direction as shown in Fig. [5] This makes the standard 
mean-field approximation fail, but may even be desirable 
[4l| . Classical chaos has a much more dramatic effect: An 
initial coherent state will rapidly spread over the entire 
chaotic sea. Thus, the description by a sharply peaked 
distribution instantaneously fails and the Bogoliubov ap- 
proach breaks down immediately. Again this breakdown 
is resolved by the introduction of a phase space ensem- 
ble as shown in Fig. QXJ The classical ensemble well de- 
scribes the spreading of the wave packet in the chaotic 
sea, demonstrating the power of the phase space ap- 
proach also for classically chaotic dynamics. 

The differences between regular and chaotic classical 
dynamics is furthermore illustrated in Fig. 1121 where we 
have plotted the dynamics of the population in the sec- 
ond well (fi2(t))/N and the magnitude of the Bloch vec- 
tor |(J(f))|/A for two different initial states. The figures 
compare the exact many-particle dynamics (solid blue 
line) to the results of an ensemble simulation (dashed 
green line) , showing a very good agreement for both regu- 
lar and chaotic dynamics. If the initial state is located on 
an island as in Fig.[T2](a) and (b), the population imbal- 
ance oscillates regularly and the magnitude of the Bloch 
vector remains close to its maximum value |(J)| = N/2. 
In this example we have chosen an initial state localized 
on a period-double fixed point (cf. Fig. ITU]) , so that the 
population difference oscillates with only half of the driv- 
ing period. The dynamics of a state initially localized in 
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FIG. 12: (Color online) Quantum dynamics in a driven 
double-well trap for an initially coherent state localized at 
(p, q) — (l/2,37r/2) on a regular island (a),(b) and localized 
at (p, q) = (0.7,0) in the chaotic sea (c),(d). The population 
in the second well {fi2(t))/N is plotted on the left-hand side 
and the magnitude of the Bloch vector | (J (t) ) | / is plotted on 
the right-hand side, both calculated exactly (solid blue lines) 
and with a classical ensemble simulation (dashed green lines). 
The remaining parameters are UN = 5, Ao = 1, to — 2tt, 
Ai = 0.5 and N = 50 particles. 



the chaotic sea is shown Fig.[T5](c) and (d). Oscillations 
are rapidly damped out and the magnitude of the Bloch 
vector drops to very small values. 

As mentioned above, the classical phase space struc- 
ture does not only affect the quantum dynamics, but also 
the organization the eigenstates (cf. Sec. IIIip . Due to the 
mixed phase space, this leads to much richer structures in 
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FIG. 13: Husimi representation of four typical Floquet states 
for a driven double-well trap with parameters UN — 5, Ao — 
1, ui = 2tt, Ai = 0.5 and N = 50 particles. 



the driven case. Instead of the eigenstates of the Hamil- 
tonian, we now have to consider the Floquet states, which 
are defined as the eigenstates of the time evolution oper- 
ator over one period T, 



U(T) =texp \-i I H(t')dt' 



(30) 



where T is the time ordering operator. Figure [T3l shows 
the Husimi representation of four typical Floquet states 
for Ai = 0.5. Three regular states are shown in the 
parts (a)-(c), corresponding to the symmetric, the self- 
trapping and the period-doubled fixed point, respectively. 
The period-doubled Floquet state shown in part (c) oscil- 
lates with a large amplitude during one driving period, in 
which it remains well localized. In contrast, the chaotic 
eigenstate are delocalized over the complete chaotic sea 
- a typical example is shown in part (d) of the figure. 



VI. BEYOND MEAN-FIELD: DEPLETION AND 
HEATING OF THE CONDENSATE 

In this section we show that the Liouville picture 
also provides an excellent quasi-classical approxima- 
tion of many-particle quantities that are not accessible 
within the common single-trajectory mean-field descrip- 
tion, since it also takes into account approximately the 
higher moments of the quantum state. In particular one 
can readily determine the deviation from a pure BEC and 
calculate the higher moments of the angular momentum 
operators ([3]) that can be used to quantify many-particle 
entanglement [4l[. In this paragraph we will especially 
discuss the first issue for the two- and three-mode case in 
more detail. This is closely related to the preceding sec- 
tions, as it is a clear signature of the regular and chaotic 
dynamics in the corresponding quantum system. 



However, the Liouville description is not restricted to 
standard Bose-Hubbard systems. Therefore, we consider 
the generalization of the two-mode system to an open 
system in the last paragraph of this section. Here, one 
quantity of particular interest is the coherence factor be- 
tween the two modes which decreases during the heating 
process [9]. It is shown that the Liouville description 
provides a valuable approximation of this decoherence 
process. 



A. The two-mode Bose-Hubbard system 

The deviation from a pure BEC can be characterized 
by the eigenvalues of the reduced single-particle density 
matrix (SPDM) 



Prcd 



1 

N 



(a\ai) (a|o 2 ) 
(ajai) 



{a\a 2 ) 



which are easily found to be given by 



2 N 



(31) 



(32) 



The leading eigenvalue of the SPDM gives the conden- 
sate fraction, i.e. the relative population of the conden- 
sate mode [121 ■ Thus, a BEC is pure if this eigenvalue 
is equal to one, i.e. if the Bloch vector lies on the sur- 
face of the Bloch sphere, |(J)|/iV = 1/2. As discussed 
above, the expectation value of the angular momentum 
operators can be calculated from the quasi-classical phase 
space distribution function p(p, q) in a good approxima- 
tion. From these it is also possible to reconstruct the 
SPDM and its eigenvalues completely classically. Fur- 
thermore, it is important to note that these quantities 
are not accessible within a single trajectory mean-field 
approach, which always assumes a pure BEC and thus 
guaranties |(s)| = 1/2 = const. To illustrate our point 
further we compare the eigenvalues of the SPDM from 
the quasi-classical Liouville approximation with the ex- 
act quantum results in Fig. 1141 
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FIG. 15: (Color online) Classical and quantum dynamics in a 
three-mode system for an initial state in the regular region of 
classical phase space (tpi(0) = 1). The upper panel shows the 
occupation in the first well Pi(t) calculated within the stan- 
dard mean-field approximation. The lower panel shows the 
quantum expectation value {hi(t))/N (solid blue line) and 
the ensemble expectation value (dashed green line). The pa- 
rameters are ei = 2, eg = 4, A12 = A23 = 1, UN — —10 and 
N = 80 particles. 
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FIG. 16: (Color online) Classical and quantum dynamics in 
a three-mode system for an initial state in the chaotic region 
of classical phase space. The upper panel shows the occu- 
pation in the third well pz(t) calculated within the standard 
mean-field approximation for pz(Q) = 1 and P3(0) = 0.999, 
respectively. The lower panel shows the quantum expectation 
value (hz(t))/N (solid blue line) and the ensemble expecta- 
tion value (dashed green line) for ^3(0) = 1. The parameters 
are ei = 2, e 3 = 4, A12 = A23 = 1, UN = -10 and TV = 80 
particles. 



Initially one eigenvalue is unity, indicating a pure BEC 
equivalent to a product state with every particle in the 
same mode. This eigenvalue decreases rapidly when the 
Husimi function approaches the hyperbolic fixed point 
(cf. Fig. [5]), indicating a rapid depletion of the conden- 
sate mode. One observes that the classical calculation 
reproduces this depletion of the condensate mode very 
well. 



B. The three- mode Bose-Hubbard system 

The same method can now be used to analyse the 
three-mode case. Two examples are shown in Figs. 1151 
and 1 161 where the standard mean-field approximation on 
the left-hand is compared to the results of a full quantum 
calculation and a classical ensemble simulation. The fig- 
ures show the population in the first and the third mode, 
respectively. Regular oscillations with a small amplitude 
are found for an initial state localized in the first well 
(V'i(O) = 1, cf. Fig. [T3]). These oscillations are weakly 
damped in the quantum system due to a dephasing pro- 
cess. This damping is well reproduced by a classical en- 
semble calculations, while there is no damping at all in 
the standard mean-field approach (upper panel). 

The dynamics is completely different for an initial state 
in the chaotic region of classical phase space as shown in 
Fig. 1161 The upper panel shows the mean-field trajecto- 



ries for two slightly different initial states (ps(0) = 1 and 
P3(0) = 0.999). They deviate fast and differ completely 
for t > 0.5. The dynamics of the corresponding quantum 
system is shown in the lower panel for an initial state 
which is completely localized in the third mode, 

|*(0)) = -^(4)^10,0,0). (33) 

The oscillations of the density (n.3) are rapidly damped 
out. Again this is well reproduced by a classical ensem- 
ble calculation. In this case the only signature of chaos 
in the common mean-field description is the sensible de- 
pendence on the initial state, however again there is no 
indication of the damping process. 

As we have already seen in the preceding paragraph, 
the ensemble simulation gives a good estimate of many- 
particle quantities, which cannot be calculated in the 
standard mean-field theory. Figure [T7] shows the evo- 
lution of the eigenvalues of the SPDM, which is defined 
analogously to Eqn. (f3"Tj) for the two-mode case. The re- 
sults from a numerical solution of the many-particle dy- 
namics are plotted as red lines, the ensemble estimates as 
dashed blue lines. The leading eigenvalues of the SPDM 
gives the condensate fraction of the many-particle quan- 
tum state. The condensate mode is depleted if the classi- 
cal dynamics is unstable in favor of the remaining modes 
[22I [23l . l26l . . The differences between the regular and 
chaotic case are obvious. As expected, the condensate 



14 



(a)1 




FIG. 17: (Color online) Evolution of the eigenvalues of the 
SPDM for (a) an initial state in the regular region of clas- 
sical phase space (V>l(0) = 1) and (b) in the chaotic sea 
(i/>3(0) = 1). Quantum results (blue solid lines) are compared 
to the classical ensemble estimate (dashed green lines). The 
parameters are chosen as in Fig. [15] and 1161 respectively. 



mode is significantly depleted if the dynamics is chaotic. 
The leading eigenvalue decreases fast, whereas it remains 
close to one if the dynamics is regular. Furthermore one 
observes a very good agreement of the classical and the 
quantum results. Therefore the ensemble method pro- 
vides a valuable tool to estimate the depletion of the con- 
densate mode and to infer the nature of the many-body 
quantum state from a purely classical calculation. This 
is especially useful for extended lattices, where quantum 
calculations become a hard problem. Another method to 
calculate the depletion of the condensate mode from the 
mean-field dynamics was introduced in [13, How- 
ever, its use is technically more involved than the simple 
ensemble approach, while it gives no additional informa- 
tion. 



C. Heating of a two-mode BEC 

Recently, the long-time dynamics of a BEC interacting 
with the b ackg round gas has attracted a lot of interest 
0, 0, l45j |. It was shown that collisions with the 
background gas lead to a decrease of the coherence of 
the two condensate modes, which can be used as a noise 
thermometer at extremely low temperatures. In this sec- 
tion we want to discuss the heating of a BEC within the 
quasi-classical phase space picture. 

The main source of decoherence and heating is caused 
by elastic collisions with the background gas atoms and 
can described in leading order by the master equation 

Ed, S3 



-i[H,p] 
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J=l,2 



n^p + ph; 



2hjphj. (34) 



The latter term describes elastic scattering events that 
conserve the number of particles in the condensate mode 
and only lead to phase decoherence (see, e.g., [11]). Thus 
they are readily described within the number conserv- 
ing phase space approach discussed in the present paper, 
where the interpretation as phase noise becomes espe- 
cially clear. 

The effects of the scattering events are conveniently 
understood and visualized in quantum phase space. The 
evolution equation for the Husimi function is obtained 
as in Sec. HIT] by taking the expectation value in SU (2) 
coherent states: 



~Q(p,g) 
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(35) 



J=l,2 

Using the I?' -algebra representation of the operators in- 
troduced above, this equation can be cast into the form 
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(36) 



7i 
" 2 



' Q(p, q) 



7=1,2 

-{H(p,q) 7 Q(p,q)} 

-2U P (i - p)q^Q(p> q) + ti^Q(p, «), 



with the classical Hamiltonian function ([23]) . By an anal- 
ogous calculation one finds the evolution equation for the 
Glauber-Sudarshan distribution: 



-{H(p,q),P(p,q)} 



(37) 



+2UP{1 ~ P) ^q~ P{P > q) + 71 W P[Pl ql 



keeping in mind that the macroscopic interaction 
strength is now given by g = U (N+2) instead of g — UN. 
In these representations the effect of the decoherence 
term ~ "fi becomes most obvious: It leads to a diffusion 
of the relative phase of the two condensates and thus to a 
blurring of the coherence of the two condensates modes. 
This effect has been directly measured in the experiments 
of the Oberthaler group [3, Eil • 

If we neglect the quantum noise term ~ g/N, the equa- 
tions ([3T]) and ([3"5|) reduce to Fokker-Planck equations. 
Thus the condensate dynamics can again be interpreted 
in terms of phase space ensembles, now subject to the 
stochastic evolution equations 



P 



dU 

-— and 
oq 



(38) 



where describes uncorrelated white noise: 

(£(*)> =0 and (Z(t)at')) = 5(t-t'). (39) 
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FIG. 18: (Color online) Heating of a two-mode BEC in quan- 
tum phase space for UN = A = 1, e = and 71 = 0.01. 
(a) Comparison of the exact quantum evolution of the Husimi 
function Q(p,q,t) (left-hand side) and the stochastic dynam- 
ics of a classical phase space ensemble (right-hand side) for 
t — 0,30,60 (from bottom to top). Dark colors encode high 
values of the Husimi density, (b) Decay of the coherence 
factor a(t). The exact quantum result (solid blue line) is 
compared to ensemble simulations based on the Q-function 
(dashed green line) and the P- function (dash-dotted red line). 



An example for this quasi-classical description of phase 
diffusion due to heating of the condensate is shown in 
Fig. [EHKa) for UN = and Fig. \W(&) for UN = 10, re- 
spectively. The left hand side shows the exact quantum 
evolution of the Husimi distribution Q(p, <?, t) calculated 
from the evolution of the density matrix p(t) according to 
the master equation At t — the system is assumed 
to be in a coherent state \p = 1/2, q = 0), i.e., a pure 
condensate with equal population and zero phase differ- 
ence between the two wells. In the course of time, the 
Husimi function spreads and thus the coherence factor 
a(t) = jj(J x )t decays as shown in part (b) of the figures. 




°' 2 20 40 t 60 80 100 
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FIG. 19: (Color online) As Fig. [18] however for UN = 10. 



The right-hand sides of Fig. UST a) and Fig. [TWa) show 
the evolution of a classical phase space ensemble initially 
distributed according to the Husimi function Q(p, q,0). 
The single trajectories evolve according to the stochas- 
tic equations of motion ([38]) and thus diffuse over the 
classical phase space. One observes that the quantum 
dynamics is well reproduced by the classical approach, 
especially the different shape of the Husimi function for 
U = and UN = 10 after the diffusion process. 

For a more quantitative analysis we have plotted the 
coherence factor a(t) as well as the fluctuations of the 
population imbalance in Fig. [T5] and Fig. 1191 The quan- 
tum result jf(J x )t is compared to the ensemble averages 
([20]) over 500 classical trajectories for the Q-function and 
P-function, respectively. In the linear case UN = the 
mapping to stochastic evolution equations is exact and 
thus the deviations only result from the finite number of 
classical representations. For UN = 10, however, quan- 
tum noise is neglected leading to a systematic underesti- 
mation of loss of phase coherence. 
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Furthermore one observes that the coherence factor 
a(t) decreases much slower for UN = 10 compared to 
the non-interacting case UN = both in the exact and 
in the classical calculation. Similarly, the fluctuations of 
the population imbalance grow much slower in the in- 
teracting case. These differences are readily understood 
from the structure of the classical phase space as shown 
in Fig. [3] One observes that the minimum of the clas- 
sical Hamiltonian function Ti.{p, q) is much deeper for 
UN = 10, so that the classical trajectories are bound 
much stronger and phase diffusion is significantly reduced 
(cf. also [3y|). In addition, this leads to a stronger re- 
striction of the phase space distribution in the J z -axis. 

VII. COMPARISON TO OTHER MEAN-FIELD 
APPROACHES 

In the preceding sections we have introduced the 
number-conserving phase space description of the Bosc- 
Hubbard model. This approach provides a straightfor- 
ward derivation of the (discrete) Gross-Pitaevskii (GP) 
equation in the macroscopic limit and allows to go be- 
yond mean-field theory. In the following we compare this 
method to other approaches to derive a mean-field ap- 
proximation. The most common ones are the Bogoliubov 
and related approximations which are discussed in detail 
in Sec. IVH Al Phase space methods are also frequently 
used, however starting from Glauber coherent states in 
almost all cases. A detailed comparison of the number 
conserving and the Glauber phase space description is 
given in Sec. IVIIBI Furthermore a considerable amount 
of work has been devoted to a rigorous derivation of the 
GP energy functional and the GP dynamics in the macro- 
scopic limit (see, e.g., (49l |) 

A. Bogoliubov theory and related approaches 

The most common method to derive the (discrete) GP 
equation is the Bogoliubov approach, taking the expec- 
tation value of the Hciscnberg equation of motion for the 
annihilation operator [50j 

ij t (a 3 ) = -J (vVi) + (a j+1 )) + U (aj.&A) (40) 
and truncating the correlation functions 

{a]a k ai) ps (ot)(ofc)(o^). (41) 

Quantum fluctuations are completely neglected in this 
equation. They can be calculated approximately with 
the Bogoliubov-dc Gennes equations which linearize the 
equations of motion for a,j around the mean-field approx- 
imation given by the GP equation [2^, [5l[ . If the initial 
state is a pure BEC, it is thus possible to calculate expec- 
tation values to all orders and to infer the increase of the 
non-condensed atoms from the dynamical (in)-stability of 



the GP equation [23J. However, this approximation in- 
cludes no backreaction of the non-condensed atoms onto 
the condensate mode, so that this approximation is re- 
stricted to short times. 

To overcome this problem, one has to truncate higher 
order expectations values at a later stage than in 
Eqn. gj), 

(ojojOj) ps (al)(a k ae) + (a]a k )(ai) + (alai)(a k ) 

-2(&})(a k )(& t ), (42) 

and to derive equations of motion for the correlation func- 
tions (ala k ) and (aja k ). Depending on the fact if and 
to which extend anomalous terms are neglected the re- 
sulting models are known under the name Hartree-Fock- 
Bogoliubov (HFB), HFB-Popov or Griffin (see, e.g., [52j 
and references therein). However, all these methods face 
some characteristic difficulties which are direct conse- 
quences of the breaking of the U(l) symmetry (53j . For 
instance, none of them conserves the total number of par- 
ticles and, at the same time, allows for particle exchange 
between the condensate and the remaining modes. 

To overcome these problems, Anglin, Vardi and Tikho- 
nenkov pro posed the Bogoliubov backreaction (BBR) 
method [26t l39l |53| . In contrast to the U[X) symmetry 
breaking approaches, only the number conserving opera- 
tors Ej k = &j«fc are taken into account. Again, equations 
of motion are calculated for the expectation values {Ej k ) 
and the correlation functions (Ej k E£ m ) , where higher or- 
der expectation values are truncated. Several numerical 
examples shown in 53] suggest that the BBR method 
provides a better approximation to the many-particle dy- 
namics than HFB and is varieties. 

But still, BBR is limited to the first two moments of the 
operators Ej k . Higher order correlation functions are not 
defined at all. This is clearly different within the phase 
space description, which allows to approximate arbitrary 
momenta. 

From a practical viewpoint, it is thus fair to say 
that the phase space method embodies the best of both 
worlds: As the Bogoliubov-de Gennes approximation, it 
allows the calculation of expectation values to all orders 
and links the increase of the non-condensed atoms to the 
dynamical stability of the 'classical' GP equation. How- 
ever, it also describes the depletion of the condensate 
mode and the long time dynamics just as the HBF and 
the BBR method. Last but not least it is even simpler 
to use and implement then both of them. 



B. Comparison to [/(l)-symmetry breaking 
methods 

Quantum phase space descriptions are usually based 
on the celebrated Glauber coherent states, which are su- 
perpositions of different number states. In particular, 
the dynamics of ultracold atoms in optical lattices has 
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FIG. 20: (Color online) Dephasing of an initially coherent 
state. Upper panel: Illustration of the dynamics of the SU(2)- 
Husimi function. Lower panel: Dynamics of the expectation 
value {L y )/N calculated exactly (solid blue line) and quasi- 
classically from an S'(7(2)-Husimi distributed phase space en- 
semble (dashed red line) and a Glauber-Husimi distributed 
ensemble (thick green line). 



been discussed in [24(. While this description is nat- 
ural in quantum optics, where the birth and death of 
photons have to be taken into account, it does not take 
into account the SU{M) symmetry of the Bose-Hubbard 
Hamiltonian. In this section, we will show that the num- 
ber conserving description is superior in describing a BEC 
with a fixed number of particles. Note, however, that the 
SU(M) symmetry is explicitly broken in the Gutzwiller 
approximation so that Glauber states provide the correct 
starting point in this case (25j . In this approximation, one 
considers only a single lattice site which couples to the 
mean-field on the remaining lattice sites so that the local 
particle number is not conserved any longer. 

The U(l) symmetry breaking phase space description 
is based on the Glauber coherent states 



\a u a 2 ) = e-CI-il J +M 2 )/2 



oo 



"1,12=0 



|m,n 2 ),(43) 



where the particle number follows a Poissonian distribu- 
tion. Let us briefly discuss the relation between the phase 
space representations based on SU (2) and Glauber coher- 
ent states. 

The Glauber coherent states (j4"5)) can be decomposed 



\ai, a 2 ) 



= P. 2 



E 

N=0 



\p> i) 



N ■ 



(44) 



into SU (2) coherent states with different particle num- 
bers N, where the respective parameters are related by 



ot\ = ae %x y/p and a 2 



pe 



(45) 



Thus, a is the amplitude and % is the global phase of the 
Glauber coherent state. Using this decomposition one 
can easily relate the Husimi functions based on SU(2) co- 
herent states and Glauber states. Consider an arbitrary 
many-particle quantum state p with N particles, which 
is represented by the SU(2) Husimi function Qiq{p ) q). 
The Husimi function in terms of Glauber states is then 
given by 



; Glauber 



(ai, a 2 ) 



[ai,a 2 \ p\oci,a 2 ) 

„-o 2 „,2JV 



JV! 



-Q N (p,q) (46) 



The physical implications of this mapping are clear: First 
of all, the Glauber-Husimi representation introduces ar- 
tificial amplitude noise, since |ai| 2 + \a 2 \ 2 = a 2 is not 
fixed but Poisson-distributed. The dynamical effects of 
this will be discussed below. The classical phase space 
is then given by C 2 and the distribution functions thus 
depend on four variables. In comparison, the approach 
based on the SU (2) coherent states has some conceptual 
advantages as it directly embodies the conservation of 
the particle number TV and eliminates the global phase 
X- While the particle number is an operator in the U(l)- 
symmetry breaking approach, it is a parameter in the 
SU{2) mapping. Thus, the discussion of the macroscopic 
limit is much simplified, as well as the estimation of the 
error, which vanishes as 0(l/N). Another approach dis- 
cussing the convergence to the mean-field approximation 
by an expansion in terms of the inverse particle number 
1/N has been used in [H,[23j]. 

The effects of the artificial amplitude fluctuations can 
be quite significant. Let us illustrate them with a sim- 
ple example shown in Fig.[20l We consider the dynamics 
of an 5(7(2) coherent state located at (p,q) = (0.7,0) 
at t = for A = 1, U = 0.1 and N = 50 parti- 
cles. This state then rotates around the elliptic fixed 
point of the classical dynamics at (p, q) = (0.5, 0) and so 
the expectation values of the angular momentum oper- 
ators oscillate. These oscillations, however, are sub- 
ject to dephasing and damping because of the quantum 
fluctuations of the initial state. Fig. [20] shows the dy- 
namics of s y = (L y )/N calculated exactly (solid blue 
line) and quasi-classically by the propagation of clas- 
sical phase space ensembles. The ensembles were dis- 
tributed according to the SU (2) Husimi function (dashed 
red line) and the Glauber-Husimi function (T4"6"|) (thick 
green line). While the damping is well reproduced by 
the SU(2) ensemble, the artificial amplitude fluctuations 
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of the Glauber-Husimi ensemble cause significant devi- 
ations. This is due to the fact that the oscillation fre- 
quency around the elliptic fixed point strongly depends 
on the effective nonlinearity Ua and thus trajectories 
with different normalization a dephase rapidly. 

Moreover, already the description of a BEC in terms 
of £/(l)-symmetry breaking states reveals some ambigui- 
ties. For instance, consider an SU (2) coherent state, rep- 
resenting a pure BEC. In the framework of the number- 
conserving phase space description, the P-representation 
of this state is simply a delta function. On the con- 
trary, the P-function in terms of Glauber coherent states 
is strongly singular, including higher-order derivatives of 
the delta-function. A pure BEC with a fixed number of 
particles is thus the most classical state within the SU (2) 
phase space description whereas it is highly non-classical 
in terms of the common Glauber phase space picture. In 
the first case an ensemble representation yields no de- 
phasing at all while such a mapping is simply impossible 
in the latter case. 

An extended and much more detailed overview over 
many different approaches to describe weakly-interacting 
condensates at finite, as well as at zero temperature is 
given in [54| . including symmetry-breaking methods, as 
well as number conserving approaches. In this terminol- 
ogy, the method presented here can be seen as a combi- 
nation of a stochastic phase space description with the 
number conserving approach. 



VIII. CONCLUSION AND OUTLOOK 

In the present paper we have discussed the number- 
conserving phase space description of small Bose- 
Hubbard systems and their generalization to an open 
quantum system. Apart from their theoretic value these 
systems have attracted considerable experimental inter- 
est within the last years, especially the two-mode case 
and its open counterpart @, @] • 

We have demonstrated the advantages of the (anti)- 
normal ordered quasi phase space densities based on the 
SU(M) coherent states instead of the common Glauber 
states for the analysis and discussion of quantum states 
and especially of their dynamics. These representations 
allow a straightforward comparison of the many-body 
quantum system and its macroscopic counterpart given 
by the celebrated Gross-Pitaevskii equation. For instance 
the quantum eigenstates localize on the classical phase 
space trajectories. Moreover, we have presented an intu- 
itive way to go beyond the usual mean-field limit: Con- 
sidering the macroscopic limit N — > oo with g — UN 
without assuming a maximally localized state yields a 
classical Liouvillian flow for the time evolution of the 
quasi-probability density or, equivalently, Monte Carlo 
ensembles of phase space trajectories. In contrast to the 
usual dynamics described by the GPE, this approach en- 
ables us to take into account arbitrary initial states and 
to estimate quantities depending on higher moments of 



the quantum state. Obviously this also holds for larger 
lattices, for which a simulation of the full many-body dy- 
namics is hard. 

As an example for the extended scope of application, 
we consider a BEC approaching a classically unstable 
fixed point. The condensate fraction rapidly decreases so 
that the dynamics cannot be described by a single GPE- 
trajectory any longer. This failure has been denoted as 
the breakdown of the mean-field approximation in the 
literature [26|, [3^|. However, this breakdown is resolved 
using the Liouville dynamics, which provides a quite rea- 
sonable approximation for higher moments. Thus, in con- 
trast to the common mean-field approach, it can also be 
used to describe and analyse chaotic systems. To under- 
line this point, we study a driven two-mode system, as 
well as the dynamics of the three-mode case, which fea- 
tures an interesting transition from a regular to a mixed 
chaotic system, depending on the interaction strength. A 
comparison to exact results shows an excellent agreement 
with the Liouville dynamics, while the common mean- 
field dynamics strongly diverge from the exact results and 
fail to describe the effects of the dynamical instability. 
Moreover, the approximation of higher moments via the 
Liouville dynamics allows to investigate many interesting 
quantities, like the single particle density matrix and the 
condensate fraction, which are not accessible within the 
common approach. To underline this point, we study the 
depletion and the heating of the condensate. 

In the last chapter, we give an overview over differ- 
ent methods, which can be used to derive mean-field 
systems and their possible extensions and relate these 
to the approach presented here. Special attention is 
payed upon the comparison to methods based on U{1)- 
symmetry breaking Glauber states and the implications 
of the conservation of the total particle number. We show 
that even for simple examples the artificial noise in the 
particle number can result in larger deviations appearing 
on a much shorter time scale compared to the approach 
presented here. Moreover, ambiguities in the description 
of macroscopic states and the set of problems consider- 
ing the macroscopic limit with TV being an operator can 
be tackled easily in the number conserving description 
where N is simply a parameter. This completes the for- 
mal comparison of different possibilities of trial states in 
[27L f55[| . Note, however, that the benefits of a description 
based on generalized coherent states are not restricted to 
dynamical problems, see e.g. the semiclassical descrip- 
tion of the ground state (56j and the excitation spectrum 
based on SU(M) coherent states [55l. 

The introduction of generalized coherent states and 
the corresponding phase space distributions opens the 
door for the use of semiclassical methods for the anal- 
ysis of the dynamics of the Bose-Hubbard model. The 
Liouville approach is of course not capable to describe 
generic quantum effects such as tunneling in quantum 
phase space and (self-)interference. These features can be 
reconstructed using semiclassical coherent state propaga- 
tors (see, e.g. [H, [5j| and references therein). In the first 
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part of the present work @ we have introduced the phase 
space description and calculated the equations of motion 
for an arbitrary number of modes. The generalization 
of the ensemble method to this case is straightforward 
and allows the approximate calculation of many-particle 
quantities such as the condensate fraction or higher mo- 
ments with small numerical efforts. Concrete applica- 
tions to larger systems, as well as an extended study of 
the three-mode system will be subject of a future paper. 

Moreover, the analysis of the mean-field dynamics has 
only recently been proven to be extremely useful for the 
understanding of the role of noise in two-mode Bosc- 
Hubbard systems which face many theoretically and ex- 
perimentally interesting features, like e.g. a stochastic 



resonance effect [4J, |45[. However, very little is known 
about the behaviour of larger dissipative systems, where 
the presented methods could provide an illustrative and 
easily implementable tool. 
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